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Abstract 

This paper proposes a new methodology to compute Value at Risk (VaR) for quan- 
tifying losses in credit portfolios. We approximate the cumulative distribution of the 
loss function by a finite combination of Haar wavelets basis functions and calculate 
the coefficients of the approximation by inverting its Laplace transform. In fact, we 
demonstrate that only a few coefficients of the approximation are needed, so VaR can 
be reached quickly. To test the methodology we consider the Vasicek one-factor port- 
folio credit loss model as our model framework. The Haar wavelets method is fast, 
accurate and robust to deal with small or concentrated portfolios, when the hypothesis 
of the Basel II formulas are violated. 

1 Introduction 

It is very important for a bank to manage the risks originated from its business activities. In 
particular, the credit risk underlying the credit portfolio is often the largest risk in a bank. 
The measured credit risk is then used to assign capital to absorb potential losses arising from 
its credit portfolio. 

The Vasicek model is the basis of the Basel II IRB approach. It is a Gaussian one factor 
model such that default events are driven by a latent common factor that is assumed to follow 
a Gaussian distribution, also called Asymptotic Single Risk Factor (ASRF) model. Under 
this model, loss only occurs when an obligor default in a fixed time horizon. If we assume 
certain homogeneity conditions, this one factor model leads to a simple analytic asymptotic 
approximation for the loss distribution and Value at Risk (VaR). This approximation works 
well for a large number of small exposures but can underestimate risks in the presence of 
exposure concentrations. 
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Concentration risks in credit portfolios arise from an unequal distribution of loans to 
single borrowers (name concentration) or different industry or regional sectors (sector or 
country concentration). Moreover, certain dependencies as, for example, direct business 
links between different borrowers, can increase the credit risk in a portfolio since the default 
of one borrower can cause the default of a dependent second borrower. This effect is called 
default contagion and is linked to both name and sector concentration. 

In credit risk management one is particularly interested in the portfolio loss distribution. 
Since the portfolio loss is usually modeled as a sum of random variables, the main task is 
to evaluate the probability density function (PDF) of such a sum. The PDF of a sum of 
random variables is equal to the convolution of the respective PDFs of the individual asset 
loss distributions. The evaluation of this convolution is a difficult problem analytically, is 
computationally very intensive and in full generality is impractical for any realistically sized 
portfolio. 

For all these reasons, several methods have been developed in the last years. The saddle 
point approximation due to [MarOlaj gives an analytical approximation of the Laplace inver- 
sion of the moment generating function (MGF). This method has been improved by |Mar06] 
based on conditional independence models. |Gla07] applies the methodology developed by 
[AbaOOj to the single-factor Merton model. First, the Bromwich integral is approximated by 
an infinite series using the trapezoidal rule and second, the convergence of the infinite series 
is accelerated by a method called Euler summation. They have shown that the cumulative 
distribution function (CDF) is comparatively accurate in small loss region, whereas the ac- 
curacy worsens in the tail region. This is because the infinite series obtained by the Euler 
summation is an alternating series, each term of which has a very large absolute value. 

Another approach to numerically invert the Laplace transform has been studied by 
[Hoo82j and |Ahn03] consisting in applying the Poisson algorithm to approximate the Bromwich 
integral by an infinite series, as in [AbaOOj and then use the quotient-difference (QD) algo- 
rithm to accelerate the slow convergence of the infinite series. We will refer to this approach 
as the Hoog algorithm. |Tak08j has applied this methodology to the multi-factor Merton 
model. The numerical examples presented show that in contrast with the Euler summation 
technique, de Hoog algorithm is quite efficient in measuring tail probabilities. 

In this paper, we present a novel methodology for computing VaR through numerically 
inverting the Laplace Transform of the CDF of the loss function once we have approximated 
it by a finite sum of Haar wavelets basis functions. This kind of functions have compact 
support and so make them useful to study local properties of the approximated function. 
Moreover, the CDF of the loss function is discontinuous, making more suitable this way of 
approximation. 

The remaining parts of the paper are organized as follows. In the next section we present 
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the one-factor Gaussian copula model and we define VaR as the risk measure used to quantify 
losses in Basel II Accord. In section three we present the basic theory about Haar wavelets 
basis system used for the approximation detailed in section four. Finally, we show with 
numerical examples the speed and accuracy of the new method in section five and section 
six is devoted to conclusions. 

2 Portfolio Loss and Value at Risk 

To represent the uncertainty about future events, we specify a probability space P) 
with sample space Q, a-algebra JF, probability measure P and with filtration (J r t)t>o satis- 
fying the usual conditions. We fix a time horizon T > 0. Usually T equals one year. 

Consider a credit portfolio consisting of iV obligors. Any obligor n can be characterized 
by three parameters: the exposure at default E n , the loss given default which without loss 
of generality we assume to be 100% and the probability of default P n , assuming that each of 
them can be estimated from empirical default data. The exposure at default of an obligor 
denotes the portion of the exposure to the obligor which is lost in case of default. Let D n 
be the default indicator of obligor n taking the following values 

_ ( 1, if obligor n is in default, 
[ 0, if obligor n is not in default, 

Let L be the portfolio loss given by: 

N 

L = L n , 

n=l 

where L n = E n ■ D n . 

To test our methodology we consider the Vasicek one-factor Gaussian copula model as 
our model framework. The Vasicek model is a one period default model, i.e., loss only 
occurs when an obligor defaults in a fixed time horizon. Based on Merton's firm- value 
model, to describe the obligor's default and its correlation structure, we assign each obligor 
a random variable called firm-value. The firm-value of obligor n is represented by a common, 
standard normally distributed factor Y component and an idiosyncratic standard normal 
noise component e n . The Y factor is the state of the world or business cycle, usually called 
systematic factor. 

V n (T) = v /p^F + a/1 - p n e n , 
where Y and e n , Wn < N are i.i.d. standard normally distributed. 

In case that p n = p for all n, the parameter p is called the common asset correlation. 
The important point is that conditional on the realization of the systematic factor Y, the 
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firm's values and the defaults are independent. From now on, we assume p n to be constant. 

Let us explain in detail the meaning of systematic and idiosyncratic risk. The first one 
can be viewed as the macro-economic conditions and affect the credit-worthiness of all oblig- 
ors simultaneously. The second one represent conditions inherent to each obligor and this is 
why they are assumed to be independent of each other. 

In the Merton model, each obligor n defaults if its firm-value falls below the threshold 
level T n defined by T n = $ _1 (P n ) where <&(#) is the standard normal cumulative distribution 
function and $ _1 (a;) denotes its inverse function. The probability of obligor ra's default 
conditional on a realization of Y = y is given by 



Pn(y) =nV n <T n \Y = y) = ^ 



Consequently, the conditional probability of default depends on the systematic factor, 
reflecting the fact that the business cycle affect the possibility of an obligor's default. 

Let F be the cumulative distribution function of L. Without loss of generality, we can 
assume X^=i E n — 1 an d consider 



F(x) = { I 
for a certain F defined in [0, 1]. 



F(x), if0<x<l, 
if x > 1, 



Let a G (0, 1) be a given confidence level, the a-quantile of the loss distribution of L in 
this context is called Value at Risk (VaR). Thus, 

l a = inf{l £l:P(I<I)>a} = inf{l G R : F(l) > a}. 

Usually the a of interest is very close to 1. This is the measure chosen in the Basel II Accord 
for the computation of capital requirement, which means a bank that manages its risks with 
Basel II must to reserve capital by an amount of l a to cover extreme losses. 



3 The Haar Basis Wavelets System 

Consider the space L 2 (IR) = {/ : f^°° \f(x)\ 2 dx < oo}. For simplicity, we can view this set 
as a set of functions f(x) which get small in magnitude fast enough as x goes to plus and 
minus infinity. 
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A general structure for wavelets in L 2 (IR) is called a Multi-resolution Analysis (MRA). 
We start with a family of closed nested subspaces 



... c v_ 2 c V-! cy cvicy 2 c ... 

in L 2 (IR) where 
and 

f{x) e \) < ^ > f{2x) e V j+1 . 

If these conditions are met, then there exists a function G Vq such that {(f)j,k}kez is an 
orthonormal basis of Vj, where 

(/) jtk (x) = 2 j/2 (j)(2 j x-k). 

In other words, the function 0, called the father function, will generate an orthonormal basis 
for each Vj subspace. 

Then we define Wj such that Vj + i — Vj © Wj. This says that Wj is the space of functions 
in Vj + i but not in Vj, and so, L 2 (IR) = J?. ®Wj. Then (see |Dau92j ) there exists a function 
ip G Wo such that {ipj,k}kez is an orthonormal basis of Wj, and {ipj t k}j,kez is a wavelets basis 
of L 2 (M), where 

V>i,fc(a;) = 2 j/2 i/;(2 : >x-k). 
The function if) is called the mother function, and the ipj^ are the wavelets functions. 

For any function / G L 2 (IR) a projection map of L 2 (M) onto V m , 

V m : L 2 (R) - V m , 

is defined by 

m— 1 fc=+oo 

V m f(x)= 22 dj,ki>j,k{ X ) = y^ C -m,fc0m,fc(^)- (1) 

i=— cxd fc=— co fcez 

where = f(x)ipj t k(x)dx are the wavelets coefficients and the c m ^ = f*™ f(x)cj) rnt k(x)dx 
are the scaling coefficients. The first part in ([I]) is a truncated wavelets series. If j were 
allowed to go to infinity, we would have the full wavelets summation. The second part in 
ffl gives an equivalent sum in terms of the scaling functions (fi m ,k- Considering higher m 
values, meaning that more terms are used, the truncated series representation of our function 
improves. 



5 



Figure 1: Scaling 



? 2,3, 



and wavelet (^2,3) functions. 



To develop our work, we have used Haar wavelets (see |Dau92] ). For these wavelets, the 
space Vj is the set of all L 2 (IR) functions which are constant on each interval of the form 



[|b ^r) for a11 integers k. Then 



and 



1, if < a; < 1, 
0, otherwise, 



1, if < x < \ 



(2' 
-1, if I < x < 1, 
0, otherwise. 

The unique thing about using wavelets as opposed to Fourier series is that the wavelets can 
be moved (by the k value), stretched or compressed (by the j value) to accurately represent 
a function local properties. Moreover, 0^ is nonzero only inside the interval [^j, ^r). 
These facts will be used later to compute the VaR without the need of knowing the whole 
distribution of the loss function. 



4 Haar wavelets approximation 

Let us mention an issue regarding the CDF F defined above. Since the loss can take only a 
finite number of discrete values (2^ at most) the PDF of the loss function is a sum of Dirac 
delta functions and then, the CDF is a discontinuous function. Moreover, the stepped form 
of the CDF makes the Haar wavelets a natural and very well-suited way of approximation. 
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4.1 Laplace Transform Inversion 



Due to the fact that F e L 2 ([0, 1]) and according to the theory of MRA, we can approximate 
F in [0, 1] by a sum of scaling functions, 

2 m -l 

F(x) m ^ C m,k<Pm,k{. X )> ( 2 ) 
fc=0 

and 

2 m -l 

F{x) = lim V c mifc mjfc (a;). 

k=0 

Recall that in the one-factor model framework, if the systematic factor Y is fixed, default 
occurs independently because the only remaining uncertainty is the idiosyncratic risk. The 
MGF conditional on Y is thus given by the product of each obligor's MGF as 

N 



M L (s; Y) = E(e~ sL | Y) = JjE [ e ~ sEnDn \ Y) 

n=l 
N 

= \{[i-Pn{y)+Pn{y)e- sE -]. 



n=l 



Notice that we are assuming non stochastic LGD. Taking the expectation value of this 
conditional MGF yields the unconditional MGF, 



M L (s) = E(e- sL ) = E(E(e~ sL | Y)) 

= E(M L (s;Y))=E 



N 



Y[[l-p n (y)+p n (y) e - sE "] 

n=l 

r " 1 2 

= / 1 f [1 -Pn(y) +Pn(y)e~ sEn ] -=e-^dy. 

JR n=1 V27T 



But if / is the probability density function of the loss function then the unconditional 
MGF is also the Laplace transform of /: 



M L (s)=E(e~ sL ) 
As we have noticed before, 

2 jv 



+oo 



e- sx f(x)dx = f(s). 



(3) 



f( x ) = ^2/J-iS(x - Xi), x 1 ,x 2 ,...,x 2 n e [0,1]. 



(4) 



i=0 



where 8{x — Xi) is the Dirac delta at Xi that can be thought as a density distribution of a 
unit of mass concentrated in the point x^ (i.e. f^°° g{x)5{x — Xi)dx = g(xi), for every test 
function g(x)). Probabilistically, a distribution, such as Q, corresponds to a situation where 
only the scenarios Xi,x 2 , ...,x 2 n are feasible with respective probabilities /ii,/U 2 , ...,{i 2 N - Of 
course these probabilities must be positive and sum up 1, this is, 

2 jv 

= i. 

As it is also well known, in the context of generalized functions, the derivative of the 
Heaviside step function is a Dirac delta. In this context (and of course in the context of 
regular functions) we can integrate by parts the expression (|3| and using the approximation 
(|2| to conclude that, 

M L (s) = / e- sx F'(x)dx = e~ s + s / e~ sx F(x)dx 
Jo Jo 

2 m -l 



■1 

e' s + s ' 



k=0 
2 m_ 1 



dx 



(5) 



e S + 2™S 22 C m,fc0m,fe(s), 



fc=0 

where ^ 

0m,fc(s) = -e" s ^(l - e~ s ^). 
s 

is the Laplace transform of the basis function (j) m ,k{x). 

— — fc 1 

Observing that = 0m,o( s ) e ~ S5,Tr an d making the change of variable z = e~ s 2^, 

the expression (JsJ) is the same as 

w^'^- ' ft( -y ) r ,r - « 6) 

2^(1-2) 

Where we note that for r < 1, Q(z) is analytic inside the disc of the complex plane {z : |#| < 
r}, since the singularity in z = is avoidable. Then, given the generating function Q(z), we 
can obtain expressions for the coefficients c mt k by means of the Cauchy's integral formula. 
This is, 

c mk = —[ ( ^-dz, k = 0,1,. ..,2 m - 1. 
where 7 is a circle about the origin of radius r, < r < 1. 



Making the change of variable z = re m , < r < 1, 

1 C 27r 0(re iu ) 1 C 27T 

Cm ' fc = 2^J ^^ du = 2^J [ Re (Q(re iu ))cos(ku) + Im{Q{re m ))sin{ku)] du 

= ~k J Re(Q(re lu ))cos(ku)du. 

(7) 

Finally, we can calculate the integral in ^ approximately by means of the trapezoidal 
rule to obtain the coefficients. 

4.2 VaR computation 

It can be easily proved that 

0<c m , fe <2-? fc = 0,l,...,2 m -l, 

and 

< C m> o < C m> i < ... < Cro^m-l. 

VaR can now be calculated quickly with at most m coefficients for each fixed level of res- 
olution m of the approximation, due to the compact support of the basis functions. Observe 
that 

— — , _ \ 22i k 1 k ~~\~ 1 

F(l a ) = 2 2 . Cm , k —<l a < — 

Thus, we can simply start searching l a computing F( 2 ^-) with any bisection like method. 

„ -=r , nm-1 » - -= , nm-1 nm — 2 » _ . — , om — 1 _i_om — 2 » _ 

If b \-™r) > a then we compute b { ^ )> otherwise we compute b ( ^ ), and 

so on. Finally, observe that VaR lies between ^ and for a certain k E {0, 1, 2 m ~ 1 }. 

5 Numerical Examples 

In this section we present a comparative study to calculate VaR between the Wavelet Ap- 
proximation (WA) method and Monte Carlo (MC). As it is well known, MC has a strong 
dependence between the size of the portfolio and the computational time. As the size in- 
creases, MC becomes a big time consuming method. 

The real situation in some financial companies show us that there are strong concen- 
trations in their credit portfolios. Basel II formulas to calculate VaR are supported under 



9 



1.E-02 T 



O.E+00 




0.15 0.20 
Loss level 



Figure 2: Tail probability approximation of portfolio PI with m — 7. 

unrealistic hypothesis, such as infinite number of obligors with small exposures. 

For these reasons, we test our methodology with concentrated portfolios. We consider 
four portfolios ranging from 100 to 10000 obligors with the main numerical results displayed 
in Table [TJ The Wavelet Approximation with m = 10 provides accurate results in a few 
seconds of computational timd^ since the maximum relative error (max|WAWA| /MC) in 
portfolios PI, P2 and P3 is only 0.5%, 1% and 0.7% respectively. The relative error in 
portfolio P4 may look somewhat greater than expected (2.5%). But in this case, performing 
5 x 10 5 MC simulations we find that the VaR obtained is 0.160, being the maximum relative 
error of 1.2% and showing again the fast convergence of the WA methodology. 



Finally, in order to display the accuracy of Wavelet Approximation, we have considered 
the tail probabilities of the loss function associated with portfolio PI using different resolu- 
tions (m = 7,8,9 and 10). The results, compared again with Monte Carlo simulations, can 
be seen in figures 2fl3f4| and [5] respectively It is remarkable how the Haar wavelets are natu- 
rally capable of detecting jumps in the cumulative distribution, making the approximation 
very precise with not many terms. 



Computations have been implemented in C under a PC with Intel CPU 280 GHz and 496 MB RAM. 
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Portfolio 


N 


Pn 


WA (m = 9) 


WA (m = 10) 


MC 


PI 


100 


0.21% 


[0.193,0.195) 


[0.194,0.195) 


0.194 


P2 


1000 


1.00% 


[0.191,0.193) 


[0.192,0.193) 


0.194 


P3 


1000 


0.30% 


[0.139,0.141) 


[0.140,0.141) 


0.141 


P4 


10000 


1.00% 


[0.160,0.162) 


[0.161,0.162) 


0.158 



Table 1: Results of 99.9% VaR computation with Wavelet Approximation and Monte Carlo 
simulations with 2 x 10 5 random scenarios and p = 0.15. In order to consider concentrated 
portfolios, in all cases we have taken E n = ^ where C is a constant such that Yln=i E n — 1. 
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Figure 3: Tail probability approximation of portfolio PI with m — 8. 

6 Conclusions 

We have presented a numerical approximation to the loss function based on Haar wavelets 
system. First of all, we approximate the discontinuous distribution of the loss function by 
a finite sum of Haar scaling functions, and then we calculate the coefficients of the approx- 
imation by inverting its Laplace transform. Due to the compact support property of Haar 
system, only a few coefficients are needed for VaR computation. 
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Figure 4: Tail probability approximation of portfolio PI with m — 9. 



We have shown the performance of the numerical approximation in four sample portfo- 
lios. These results among other simulations, show that the method is applicable to different 
sized portfolios and very accurate with short time computations. Moreover, the Wavelet 
Approximation is robust since the method is very stable under changes in the parameters 
of the model. The stepped form of the approximated distribution makes the Haar wavelets 
natural and very suitable for the approximation. 

We also remark that the algorithm is valid for continuous cumulative distribution func- 
tions, and that it can be used in other financial models without making conceptual changes 
in the development. For instance, we can easily introduce stochastic loss given default (just 
changing a bit the unconditional moment generating function) and consider the multi-factor 
Merton model as the model framework as well. 
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